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Abstract 

A pulsating form of hydrodynamic instability 
has recently been shown to arise during the deflagra- 
tion of liquid propellants in those parameter regimes 
where the pressure-dependent burning rate is char- 
acterized by a negative pressure sensitivity. This 
type of instability can coexist with the classical cel- 
lular, or Landau, form of hydrodynamic instability, 
with the occurrence of either dependent on whether 
the pressure sensitivity is sufficiently large or small 
in magnitude. For the inviscid problem, it has been 
shown that when the burning rate is realistically al- 
lowed to depend on temperature as well as pres- 
sure, that sufficiently large values of the tempera- 
ture sensitivity relative to the pressure sensitivity 
causes the pulsating form of hydrodynamic insta- 
bility to become dominant. In that regime, steady, 
planar burning becomes intrinsically unstable to pul- 
sating disturbances whose wavenumbers are suffi- 
ciently small. In the present work, this analysis is 
extended to the fully viscous case, where it is shown 
that although viscosity is stabilizing for intermedi- 
ate and larger wavenumber perturbations, the in- 
trinsic pulsating instability for small wavenumbers 
remains. Under these conditions, liquid-propellant 
combustion is predicted to be characterized by large 
unsteady cells along the liquid/ gas interface. 

Introduction 

Hydrodynamic (Landau) instability in combus- 
tion is typically associated with the onset of wrin- 
kling of a flame surface, corresponding to the for- 
mation of steady cellular structures as the stability 
threshold is crossed. This type of instability was 
originally described by Landau 1 , and is attributed 
to thermal expansion across a combustion front. Al- 
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though gaseous combustion was determined to be 
intrinsically unstable in Landau’s analysis, it was 
demonstrated that additional effects, such as gravity 
and surface tension, that enter when the unburned 
mixture is a liquid result in a specific stability cri- 
terion. However, this analysis, along with a later 
study by Levich 2 that considered viscous effects in 
lieu of surface tension, assumed that the the combus- 
tion front propagated normal to itself with constant 
speed, whereas it is now recognized that there is a 
dynamic interaction between the burning rate and 
local conditions at the front. 

For those problems in which pyrolysis, exother- 
mic decomposition and/or combustion occurs in an 
intrusive region in the vicinity of the liquid/gas in- 
terface, the dynamical coupling of the burning rate 
with the underlying hydrodynamics of the flow can 
be achieved through an analysis of the thin com- 
bustion/interface region. An alternative approach, 
however, is to simply postulate a phenomenologi- 
cal dependence of the local burning rate on pressure 
and temperature, and to obtain results in terms of 
suitably defined sensitivity parameters. Both types 
of methodologies have been applied to the problem 
of solid-propellant combustion, and each offers cer- 
tain advantages. 3,4 In the present series of studies 
on liquid- propellant combustion, 5 7 the latter ap- 
proach has been adopted, thereby generalizing the 
Landau/Levich model to allow for a coupling of the 
burning rate with the local pressure and tempera- 
ture fields. 

Summarizing some of the results obtained from 
the present model, it has been shown that when only 
the pressure sensitivity of the burning rate is taken 
into account, an appropriately generalized stability 
criterion for cellular (Landau) instability is obtained. 
Exploiting the realistic limit of small gas-to-liquid 
density ratios, it is found that the stable region oc- 
curs for negative values of the pressure-sensitivity 
parameter, with the original Landau model being 
intrinsically unstable in this limit. In particular, 
the neutral stability boundary possesses a local min- 
imum when it is plotted against the wavenumber 
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of the disturbance, suggesting that as the pressure- 
sensitivity parameter decreases in magnitude, the 
liquid/gas inter face/front develops cells correspond- 
ing to classical hydrodynamic instability. This min- 
imum reflects the fact that surface tension and vis- 
cosity are stabilizing influences for short-wave dis- 
turbances, whereas gravity is a stabilizing influence 
for long- wave perturbations. As a result, the effect of 
reducing the gravitational acceleration to micrograv- 
ity levels is to shift the neutral stability minimum 
to smaller wavenumbers. Thus, in the micrograv- 
ity regime, Landau instability becomes a long-wave 
instability phenomenon, implying the appearance of 
large cells along the combustion interface. 

Aside from the classical cellular form of hydro- 
dynamic instability, this dynamic generalization of 
the Landau/Levich model also predicts the appear- 
ance of a pulsating form of hydrodynamic instabil- 
ity, corresponding to the onset of temporal oscilla- 
tions in the location of the liquid/ gas interface. This 
form of hydrodynamic instability occurs for negative 
values of the pressure-sensitivity parameter that are 
sufficiently large in magnitude. 6 Consequently, sta- 
ble, planar burning is predicted to occur in a range 
of negative pressure sensitivities that lies below the 
cellular boundary and above the pulsating boundary 
just described. A stable range of negative pressure 
sensitivities is applicable, for example, to certain 
types of hydroxylammonium nitrate (HAN)-based 
liquid propellants at low pressures for which non- 
steady modes of combustion have been observed. 8 
The appearance of both pulsating and cellular forms 
of hydrodynamic instability is analogous to the two 
corresponding types of thermal/diffusive instabili- 
ties that occur for sufficiently large and sufficiently 
small Lewis numbers, respectively. 9 

When the additional effect of a temperature 
sensitivity in the burning rate is included in the 
analysis, substantial modifications to the above sta- 
bility description can occur. Specifically, it is found 
that if the temperature-sensitivity parameter is suffi- 
ciently large relative to the parameter corresponding 
to pressure sensitivity, the pulsating hydrodynamic 
stability boundary can develop a turning point (i.e., 
become C-shaped) in the (disturbance- wavenumber, 
pressure-sensitivity) plane. In that case, the sta- 
ble region for small wavenumbers disappears, and 
liquid-propellant combustion is predicted to be in- 
trinsically unstable to the nonsteady form of hydro- 
dynamic instability for all sufficiently large distur- 
bance wavelengths. This has been described in de- 
tail in the limit of zero viscosity, 7 and the purpose 
of the present work is to extend that analysis to the 
fully viscous model. Viscous effects were shown to 


have a substantial influence in the absence of ther- 
mal sensitivity, where it turned out that the stable 
region became significantly widened when viscosity 
was present, and the same result will be demon- 
strated when thermal effects are present. However, 
the same intrinsic pulsating instability that occurs 
for sufficiently large temperature sensitivities and 
sufficiently small wavenumbers in the inviscid case 
will be shown to be preserved even when viscosity is 
included. These results lend further support to the 
notion that a likely form of hydrodynamic instabil- 
ity in liquid-propellant combustion is of a nonsteady, 
long-wave nature, distinct from the steady, cellular 
form originally predicted by Landau. 

Summary of the Mathematical Model 

The present model of liquid-propellant combus- 
tion was described previously, 8 ** 6 but is briefly sum- 
marized here for completeness. Specifically, it is as- 
sumed that the combustion front coincides with the 
liquid/gas interface, where pyrolysis and/or exother- 
mic decomposition occurs. Denoting the nondimen- 
sional location of this downward-propagating inter- 
face by £3 — where 23 is the vertical co- 

ordinate and the adopted coordinate system is fixed 
with respect to the stationary liquid at 23 = -00, we 
transform to the moving coordinate system 2 = 21, 
y = x 2 , 2 = X3-${xx } x 2l t) such that the liquid/gas 
interface always lies at z = 0. Conservation of mass, 
energy and momentum within each phase then gives 

V * v = 0 , z/0, (1) 


50 

dt 

dv 

~dt 


^ +V .V 0 =(J)V 2 0 , z<0, ( 2 ) 

dt dz '■AJ 

^ + (v .V)v = (0,0,-*Y-‘> 

2 5 °' 


where v, © and p denote velocity, temperature and 
pressure, respectively, Pvi and Prg denote the liquid 
and gas-phase Prandtl numbers, p, A and c (used be- 
low) are the gas-to-liquid density, thermal diffusiv- 
ity and heat-capacity ratios, and Fr~ l is the inverse 
Froude number (gravitational acceleration). Other 
nondimensional parameters that appear below in- 
clude the liquid surface-tension coefficient 7, the gas- 
to-liquid viscosity ratio p ( p\Pr g — pPr{), the rate- 
of-strain tensor e and the unburned-to-bumed tem- 
perature ratio <r u . 
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Equations (1) - (3) are subject to a set of bound- 
ary and interface conditions given by 


v = 0 , 6 = 0 at z — -oo , 

0 = 1 at z = +oo , ©| 2 =o- = ©U=o+ i 


( 4 ) 


n s x v_ = n s x v+ , (5) 

0$ 

n,-(v_ - pv+) = (1 - p)S{$)-^ , (6) 

n s v_ — = -4(p| 2 =o+ i ©U=o) > CO 

p|z=o- p|z=o + 

= n s - [pv+( n s • v + ) - v_(n, • v_) 

-p\Pr g e+ ■ n s + Pne. • n s ] (®) 

5$ 

+ n g (v_ - pv+)S(Q)-jg - 7 (-V • n.) , 
n s x |pv + (n s • v+) - v_(n s • v_) 

+ (v_ -pv + )S($)^] 

= n s x(p\Pr g e+ • n s - Pne_ * n 8 ) , 

n 9 • (cpAV0| z= o+ - V0| z -o-) 

= n s * [(cpv+ - v_)0| 2=o + c(a u pv+ - v-)] ( 10 ) 

d$ 

+ [(1 “ cp)©| 2 — o c(l ~ (Tup)] *S(^) q- , 


where c = c/(l - a u ), v± = v| 2== o±> e ± = e lz=o ± i 
and Eqs. (5) - (10) correspond to continuity of 
the transverse velocity components (no-slip), con- 
servation of (normal) mass flux, the mass burning 
rate (pyrolysis) law, conservation of flux of the nor- 
mal and transverse components of momentum, and 
conservation of heat flux. Here, the surface fac- 
tor S($) = (1 + $1 + $J) -1/2 . the unit normal 
n 3 = ( — <J> X , -$ y , l)S($), and the expressions for 
the gradient operator, the Laplacian and the curva- 
ture in the moving coordinate system are given by 
V = (d x - $ x d 2 , 9 y - $ v d z ,d z ), V 2 = d xx + d yy + 
(1 +S* 2 + $ 2 ) d zl - 2$ x d xz - 2$ y d yz - ($ IX + $ yy )d 2 
and — V • n a = **.(l + *£) + ®» y (l + *2) - 
2$x$y^xy)i respectively. However, the vector v still 
denotes the velocity with respect to the (x,, x 2 ,x 3 ) 
coordinate system. Finally, we note that the nondi- 
mensional mass burning rate appearing in Eq. (7) is 
assumed to be functionally dependent on both the 
local pressure and temperature at the liquid/ gas in- 
terface. By definition, .4 = 1 for the case of steady, 
planar burning, but local perturbations in pressure 
and/or temperature result in corresponding pertur- 
bations in the local mass burning rate. 


Since the thermal and hydrodynamic fields are 
coupled only through the temperature dependence 
of the mass burning rate A appearing in Eq. (7), 
the strictly hydrodynamic problem for p, v and 
can be analyzed separately when A is assumed to de- 
pend on pressure only [5,6], In the present work, we 
focus on the fully coupled problem to determine how 
the hydrodynamic stability boundaries are modified 
when the local burning rate depends on tempera- 
ture as well as pressure. Our stability results will 
thus depend on two sensitivity parameters, A p and 
.4©, defined as A p — ch4/9p| e=lp=0 and Aq = 
5A/ae| 0=1 =0 , where © = 1 and p = 0 denote the 
interface values of temperature and pressure of the 
basic solution in Eq. (11) below. Though an explicit 
expression for the reaction rate A is not needed in 
the present analysis, we note that since the nondi- 
mensional activation energy is typically large, the 
temperature sensitivity Aq would likely be larger in 
magnitude than the pressure sensitivity A p , a fact 
that will have some bearing on the relative scalings 
of these parameters that will emerge in the following 
analysis. 

The Basic Solution and its Linear Stability 

The nontrivial basic solution of the above prob- 
lem that corresponds to the special case of a steady, 
planar deflagration is given by 
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>°(z) = 

1 -Fr~ 1 z + p~ x 
X-pFr-h, 

“ 1, * 
Z 

<0 
> 0 

«°), 

„o _ f 0, 

V \ P -1 “ 1. 

z < 0 
z > 0 , 


(e z - 

z < 0 



U 

z > 0. 




( 11 ) 


The problem governing its linear stability may be 
formulated, prior to introducing any further approx- 
imations, in a standard fashion by introducing per- 
turbation variables <f>(x,y,t) = $(x, y, z, t) - $°(t), 
u (x,y,z,t) = v(x,y,z,t) — v°(z), C(x,i/,z,f) = 
p{x,y,z,t) — p°(z) and 6(x,y,z,t) = 0 — ©°(z) — 
<j> a d&°/dz. The stability problem obtained when 
Eqs. (1) - (10) are linearized about the basic so- 
lution (11) is then given in terms of these variables 


by 


dui C)U2 du3 _ Q 

dx dy dz 


z/0, 


( 12 ) 
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r 1 1 9u 9u __ 

1 pi'dt + ~dz ~ 

fd<: jii _i^ 9C , rM F , -i££ 

\ax + lp/ Fr dx’dj/ + ip) 5y’9zj 

f Pr t wd 2 u a 2 u jggA 2<0 
I fiPr\ J \ dx 2 dy 2 dz 2 / ’ > 

f 1 1 90 90 _ r -U 3 e z | 

1 p J 9 £ 9z l 0 J 

r 1 \/9 2 0 ^d 2 9 9 2 0\ < 

+ {^/fe + 9? + ^’ 

(15) 


_iu>t+tfcix+tfc 2 y 


(13) 

(14) 


u = 0 , 0 = 0 at z~ -oo , 

0 = 0 at z — “boo , 0| z =o+ ** 0|z= o- = 

Ul| z= 0- -Ul| z=0 + = (^ _1 _1 )^> 
W2| z=0 - - u 2| z=0+ = (P _1 “ 1)^V > 

U3| z=0 - -PW 3 | Z=0+ = (1 “ P)fa , 
u 3 | z=0 _ - 4>t = >lpCU=o+ + ^©^U=o+ 1 

<U=o- - CU=o+ = 2 W z =o+ - u 3 | z=0 -) 


2 Pn ( 


dU3 

dz 


du3 
z =o- ^ dz 


E=0+) 


- 2(1 - p)4>t ~ 7 (<t>xx + <f> V y) . 

_ fdu\ du 3 \ \ n 

-Pni-z 1 + -s— „ ) =°> 

\ dz z= o- dx lz=o-/ 

^(SU + $U 

<9u2| , 9u3 


cpA 


-Pn(P I + ^r n )=°> 

V 9z 1 2=0- 9y 2-0-/ 

901 901 . m 

— - - c0| z=o + + V\z=o- 

9zlz=o+ oz \z=Q~ 

= -u 3 | z=0 - +0t, 


(16) 

(17) 

(18) 

(19) 

( 20 ) 
( 21 ) 


where Eqs. (16) and (17) have been used to simplify 
Eqs. (18) - (21). 

Nontrivial harmonic solutions for </>, u and C > 
proportional to e iut+iklX+ik3V , that satisfy Eqs. (12) 
- (14) and the boundary/boundedness conditions at 
z = ±oo are given by 4> = e tult+lklX+tk ' il1 and 


f _ iut+ikix+ikiy / b\e kz Fr , Z < 0 
C “ \b 2 e~ kz - pFr~\ z> 0, 


( 22 ) 


__ Au)t+ikix+ik 2 y 

u i = e 


f b 3 e qz — ifci(io; + k)~ 1 bie kz , z< 0 ( 23 ) 

\b 4 e rz -ik 1 (iujp-k)- 1 b 2 e- kz , z> 0 , 


u 2 = e 

fhe , < z -ik 2 {iu} + k)- 1 bie kz , z< 0 ( 24 ) 

' \fe 6 e rz -tA: 2 (iu;/j-fc)- 1 6 2 e- fc fy z>0, 

_ -ivt + ikix+ikiy 

U$ — e 

j b 7 e qz -k{iw + k)~ 1 bie kz , z< 0 ( 25 ) 

[b s e rz +k(icjp-k)~ l b 2 e~ kz , z> 0, 

0 __ giut+ifcix-^ifcay 

f 6ge pz - [iu + Jfc 2 - q(< ? +l)]- 1 6 7 e^ +1 ) z 

+ fc[(iw) 2 -fc 2 ]- 1 6 1 e( fc+1 )fy z<0 

* > (*> 

Here, k = (k 2 + fc 2 ) 1 ^ 2 is the overall disturbance 
wavenumber and the quantities p, Q , r, s are de- 
fined respectfully as 2p = 1 + [l + 4 (iui + A: 2 )] > 

2 Pn q = 1 + [l + 4Pn(iu + Pn A: 2 )] 1/2 , 2pPr t r = 

1 - [l + ipPri(iup + pPn A; 2 )] and 2pAs = 1 - 
[l + 4p 2 A(tw + AJt 2 )] 1/2 . Substituting this solution 
into the interface conditions (16) - (21) and using 
Eq. (12) for z > 0 yields eleven conditions for the 
ten coefficients 6i - feio and the complex frequency 
(dispersion relation) iu(k). In particular, these con- 
ditions are given by 

iAqb 3 4- ife 2 6s + qbz = iAqt^ 4- iA: 2 ?>6 + r b& — 0 1 ( 2 ^) 

b 3 - -r—^rbi - 64 + . — r^2 = (p 1 “ l) ik i < 

iuj H - k xu)p fc 

65 - . r h - be 4- 7 ^2 = (p 1 ” 1)^2 » 

to; + fc lup-k 

(28a, 6) 

^ - pb 8 - 1 — — r^2 = (1 - 0)^ 1 

iu -b /c - k 

b 7 — t&i — -Ap 62 ^10 = it*; — pF r 1 A p , 

- + * (29a, 6) 


10/ 


fl + H^£2l6, - [l + + 1 - p)]^ 

L iu; + A: J L iup-k J 

- 2Pr t qb 7 -2(1 - p - pPri r)b s 

= (1 - p)(Fr _1 - 2iu;) + 7A: 2 , 

pPr, r6 4 + 2ik ^ P p b 2 + ihpPn b s 
^ tup — k 

- Pri qb 3 H : — —7—61 - ikiPn 67 = 0 , 

to; -1- A: 

pPri rb 6 + + ik 2 pPri b a 


(30) 


(31) 


- Pri qb 5 + 


iujp — k 
2ik2kPri 


(32) 


XU) -b k 


b x - ik2Pri b 7 = 0 , 


4 


6io — ^7 A;[(tw) k ] b\ 1, 

(33) 


(1 - c + cp\s)bio + 


<? + 1 

iw + k 2 - q{q + 1) 



b 7 


-pb 9 - 


k 

iu> + k 


k + 1 
iu - k 



b\ = 1 4- iw . 

(34) 


While the above problem is linear in the coef- 
ficients b\ — 6io, which can thus be eliminated to 
give a single equation for iw, the resulting disper- 
sion relation is quite long and highly nonlinear. Ex- 
plicit results may be obtained for certain special 
cases, including the original problems considered by 
Landau 1 (A p = Aq = Pn = p = 0) and Levich 2 
(A p = Aq = p = 7 = 0), as well as a particular 
case (A© = p = Pr t = 0) of the generalized model 
described above. 5,11 To obtain more general results, 
it is possible to exploit the smallness of certain pa- 
rameters and to seek asymptotic solutions for the 
neutral stability boundaries. In particular, realistic 
limits to exploit include the smallness of the gas- 
to-liquid density and viscosity ratios p and p, and, 
in the microgravity regime, Fr -1 . Pursuing this ap- 
proach, tractable asymptotic results have so far been 
obtained for the case 5,6 Aq = 0 and for the invis- 
cid problem when >1© is nonzero. 7 The present work 
essentially completes the asymptotic analysis of the 
dispersion relation embodied in Eqs. (27) — (34) by 
extending the last of these studies to the fully vis- 
cous problem. 


Parameter Scalings and Asymptotic 
Analysis of the Dispersion Relation 


instability using the generalized model in the limit 
A e = 0, but they are also relevant when one consid- 
ers the pulsating form of instability and when Aq 
is allowed to be nonzero.In particular, in the case 
of cellular instability and zero temperature sensitiv- 
ity, there are three wavenumber scales to be con- 
sidered. First, there is an 0(1) outer wavenumber 
region where the stabilizing effects of surface ten- 
sion, viscosity and gravity are all relatively weak. 
Second, there is a far outer scale k ~ fc//e where 
surface tension and/or viscosity are strongly stabi- 
lizing and gravitational effects are, to a first approx- 
imation, negligible. Finally, we have an inner scale 
k ~ ki€ or k — M 2 where gravity is the dominant 
stabilizing effect (the first scale is valid for normal 
gravity, the latter for the reduced gravity regime), 
and where viscosity and surface-tension effects are 
absent at leading order. In each of these regimes, 
the cellular stability boundary, obtained by seeking 
solutions of the dispersion relation for which iu> is 
identically zero, is given respectively by 


A;(k)~-p‘/2, 



p*(p* 9 ~ ki)/2ki, 
P*{P* 9* ~ ki)/2k iy 


Fr -1 — 0(1) 
Fr" 1 -O(e), 

(35) 


a; u) ~ -p* 

2pTj£P_ [1 + fc/(p*7 + 2 ^* P + 2 P*F)] 

+ 4/r*P(l + p*Pfc/) - [1 - R{kf)](p* 7 + 2 M*P) ’ 

(36) 

where P = Frj and R(kf) = (1 + 4 p* P 2 k 2 ) 
Matching these solutions to one another then leads 
to the composite stability boundary 


Focusing on the realistic regime p«l (typical 
values are on the order of 10“ 3 or 10 4 ), we formally 
introduce a bookkeeping parameter e«l and intro- 
duce the reasonable scalings p = p*e y p = p* t, Pri ~ 
0(1) and either Fr' 1 = g or Fr' 1 = £r*e, where 
Fr' 1 ~ O(e) corresponds to the case of greatly re- 
duced gravity. In this parameter regime, the appro- 
priate scaling for A v to describe the neutral stability 
region is A p = A*e [5,6], whereas the appropriate 
scale that describes the main effects of thermal cou- 
pling turns out to be Aq = AqC 1 / 4 [7]. Based on this 
scaling, we note that the ratio Aq/A p ~ 0(e' 3 ^ 4 ) 

1, as might be expected based on an overall Arrhe- 
nius reaction-rate dependence on temperature. 

Based on our previous analyses, the above scal- 
ings induce a set of corresponding regimes for the 
wavenumber k (and the complex frequency iw) in 
the dispersion relation determined by Eqs. (27) - 
(34). These first emerged in our analysis of cellular 


A* (c) 


2 p'p'P [1 + £fc(p*7 + 2M*F + 2p*P)] 


tg 

eV 


+ 4/i*P(l + tkp*P) - [1 - R{ek)]{p * 7 + 2 p*P) ’ 

(37) 


as shown in Fig. 1. Clearly, the stable region lies 
below A* = —p* 1 2 (negative values of A p over cer- 
tain pressure ranges are characteristic of a number of 
HAN-based liquid propellants [8], with the location 
of the minimum in the cellular boundary increasing 
to less negative values of A p with increasing values 
of the stabilizing parameters 7*, Pn , p* and g (or 
g*). Comparing the two sets of curves corresponding 
to the normal and reduced gravity cases, it is clear 
that the critical wavenumber for instability becomes 
small in the latter regime. That is, cellular hydrody- 
namic instability becomes a long-wave instability in 
the limit of small gravitational acceleration. Further 
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discussion of this stability boundary, and its rela- 
tionship to the original Landau/Levich predictions, 
is given in Ref. [5]. 



Figure 1. Asymptotic representation of the cellular 
hydrodynamic neutral stability boundaries. The up- 
per and lower sets of curves correspond to the nor- 
mal and reduced- gravity regimes , respectively (curves 
drawn for the case e = .04, p* = 1.0, g = 6.0, 
g* = 2.0 ). The solid curves correspond to the invis- 
cid limit (P = 0) with nonzero surface tension = 
2.5j. The dash-dot curves correspond to nonzero 
surface tension (*y — 2.5) and liquid viscosity (P — 
l.o), but zero gas-phase viscosity (p*P = 0J. The 
dash-dot-dot curves differ from the dash-dot curves 
by the addition of gas-phase viscosity (fX*P = l.Q), 
and are similar to the dash-dot-dot-dot curves , where 
the latter correspond to larger viscosities (P = p* P 
= 2.0). The dash-dot-dot-dot-dot curves correspond 
to a viscous case (P = fi* P = 1-0^, but with zero 
surface tension. 



Figure 2. Asymptotic representation of the pulsating 
hydrodynamic neutral stability boundary for the vis- 
cous case (P > 0). The region between the pulsating 
and cellular boundaries (the latter are shown on an 
expanded scale in Fig . 1) is the stable region with 
respect to hydrodynamic instability . 


Considering the pulsating stability boundary in 
the limit Aq = 0, which is obtained by seeking so- 
lutions of the dispersion relation for which only the 
real part of iu vanishes, it was determined 6 that 
the corresponding expressions in the inner and outer 
wavenumber regions are given by 

A* ~ -p* , A; ~ -p*(l + 2 Pk) 1/2 , (38) 

respectively. In this case, it is clear that the outer 
solution is, in fact, the composite solution, which lies 
below the cellular boundaries and recedes to nega- 
tive values of A p that are larger in magnitude than 
O(e) as k becomes large (Fig. 2). Clearly, this sta- 
bility boundary is more sensitive to the stabilizing 
effects of the liquid viscosity parameter P than is 
the cellular boundary, having a leading-order sta- 
bilizing effect for 0(1) wavenumber disturbances in 
this case. In the limit P -* 0, the pulsating bound- 
ary collapses to the straight line A p = —p (i. e., 
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A* = -1 in Figs. 1 and 2), but even in that limit, 
there is a region of stability corresponding to val- 
ues of Ap greater than — p m and less than the mini- 
mum in the cellular boundary, which is greater than 
-p* / 2. However, if one now considers the effects of a 
nonzero temperature sensitivity in the inviscid limit 
P = 0, then, for Ae ~ C^e 1 / 4 ), the pulsating bound- 
ary possesses a turning point such that the stability 
region disappears for sufficiently small wavenumber 
perturbations. 7 This is illustrated in Fig. 3, which 
indicates that the pulsating boundary then frames 
the stable region except along the upper branch that 
asymptotes to the previous cellular boundary as k 
becomes large in the outer wavenumber region. The 
evolution from a stability diagram that indicates a 
stable region delineated by distinct pulsating and 
cellular hydrodynamic stability boundaries to the 
pulsating-dominated one exhibited in Fig. 3 can 
be shown to occur in the parameter regime Ae ~ 
0(e 1 / 2 ), which, based on the estimate Ae/A p ~ 
Oe -1 / 2 >30 (t.e., of the same order as a typical 
nondimensional activation energy), appears to be at- 
tainable for many types of liquid propellants. We 
now extend the analysis that produced the fully- 
developed pulsating boundary shown in Fig. 3 to 
the viscous case in which both P and p* are allowed 
to be nonzero. 


Owing to the complexity of the fully viscous 
problem, we analyze Eqs. (27) - (34) directly by 
seeking appropriate expansions for the complex fre- 
quency tut and the coefficients 6,. This differs from 
our approach in the inviscid limit where it was feasi- 
ble to first eliminate the bi in order to obtain a single 
implicit equation for tw alone. We first consider the 
0(1) wavenumber region and, based on our previ- 
ous analyses, seek an expansion for the dispersion 
relation tu>(fc) in this region in the form 

iw ~ c -1/,2 (tu; o + € X ^ 4 wvi + H ) i (39) 

Introducing the previously defined parameter scal- 
ings, the quantities p, q, r and s defined below Eq. 
(26) are expanded as 

p ~ p(_i/4)C -1 / 4 + Po + P(i/4) e ^ ) 

q ~ 9(- 1/4 )£ _1/4 + 9o + • • ' . ( 40 ) 

r ~ r(i/ 2 )e 1/2 + r 3 /4« + r i c d • 

s ~ S(i/ 2 )(e 1/2 ) 4 > 

where p(_i/ 4 ) = (iw o) x ^ 2 t 2po = [l + iwi/(iwo) ^ ], 
8 p (1/4 ) = (iw 0 ) -1/2 [l + 4fc2 + 4iu;2 - (iwi) 2 /iWo]. 


Composite Pulsating/Cellular Stability Boundary 


Aq = .5, p' = 7 = 1 . *=005 



Figure 3. The composite pulsating /cellular hydrody- 
namic stability boundary for Ae ~ 0(e 1 ^ 4 ) in the 
limit of zero viscosity. 
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g ( _i/ 4 ) = (wo/P) 1 ' 2 , 2 Pq 0 = [1 + tw 1 /(iwo/P) 1/2 ], 

r (i/2) = s (i/2) = — vjJqP* , *"(3/4) — p*, and r i = 

_ iu) 2 p * - ( p*Pk ) 2 . Finally, the 6, are conservatively 
postulated to have the expansions 


6, = 6|- V 1 + fc'- 3/ V 3 / 4 + b^r 1 ' 2 + • • ■ , 

• = 1,2,8, 

6 i = 6<- 1/2) £- 1 / 2 + fci - 1 / V 1 / 4 + 6 [ 0) + -" 1 

< = 3,4, 5, 6, 


6 . = 6 (-i/V 1 / 4 + 6|V + 6i 1/4) £ 1/4 + ---, 

i = 7, 9, 10, 

(41) 

where the form of the latter expansions is again 
partly motivated by our previous analyses of more 
specialized cases. 

Substituting these expansions into Eqs. (27) - 
(34) and equating coefficients of like powers of e, we 
obtain the leading-order conditions 


+M-' m +nm/r l) =<>• 


4-» _-*•//, . . 

l(-i) i i(-l) _ o// -1 ) — 0 

®1 + °2 20 8 - U - (44) 

bl'V = -(k/p')(i+ a;/ p'), 

where the last of these equations was obtained from 
the leading-order difference of the first and second 
of Eqs. (29) using the last of Eqs. (43). From Eqs. 
(43) and (44) we thus conclude that the leading- 
order dispersion relation is given by 

(iw Q ) 2 = (k/p*) 2 (2A; + p-), (45) 

which is the same result as that obtained in the in- 
viscid case when Aq = 0. In particular, Eq. (45) 
implies that (iw 0 ) 2 % 0 for A\ \ - p*/2, which re- 
covers the leading-order cellular boundary (35) for 
0(1) wavenumbers, but gives no definitive informa- 
tion regarding stability for A* < —p */ 2 because two 
is purely imaginary in that region. That is, the sta- 
bility of the basic solution in the latter region is 
determined by the real parts of higher-order coef- 
ficients in the expansion (39) for iu>, although the 
fact that Im{iu>o} £ 0 implies that disturbances 
have a pulsating character for values of A* below 
the cellular stability boundary. 


At the next order in the analysis of Eqs. (27) - 
(34), we obtain a second set of conditions given by 

iki &3 -1/4) +**2 *4 -1/4> +9(-i/4) 4 0> +9o 4 ^ 4) 


(47) 


(48) 


ihb[' l/4) +ik 2 b { 6 l/4) 

+ *■(1/2) *4 3/ ^ + r (3/4)&8 

i Wl =6i- 3/4 > = 6<- 3/4) =^ 1/2) =^ l/2) 

= 6<- 1/4) =6' O) =6 ( 1 ° o , = 6'- 3/4) =0, 

where the last of Eqs. (48) was deduced from the 
next-order difference of Eqs. (29). Finally, from the 
sum of the first of Eqs. (28) multiplied by iki and 
the second of Eqs. (28) multiplied by ik 2 , we con- 
clude that 6^" 1/2) = iw 0 (l - AJ/p*). However, the 
fact that tu>i = 0 implies the need to continue the 
analysis at the next order to determine tw 2 . Pro- 
ceeding in this fashion, we obtain from the previous 
results and Eqs. (29) - (34) at this next higher order 




rtf r'rtnrlifi/'kne 




b^-(k/iw 0 )b[- 1/2) -2iu; 2 = k(l-A;/p*), 4 0) = 0, 

(49) 

b'- 1/2) - 2<4" 1/2) = iu, 0 (2fcP - 1 + A;/p*) , (50) 


6'- 1/2) -(A^/p*)6^ 0 /4) = ^o[l-(A;/p*) 2 ] , (51) 


- (k/iu> 0 )b[ 1/2) 


- (ia> 0 ) 1/2 6io /4) 


- 2iw 2 = 2Pk 2 , 


(52) 


where Eq. (51) was actually obtained from the next 
higher order difference of Eqs. (29), and the second 
of Eqs. (49) was obtained from the sum of Eqs. (31) 
multiplied by iki and Eqs. (32) multiplied by ik 2 . 
Equations (49) - (52) constitute a closed system for 
b[~ l/2 \ 6g~ 1/2) , 6 jo / 4) and iw 2 . Eliminating the first 
three of these in favor of the last and using the result 
(45) for iw 0l the dispersion relation for iw 2 is finally 
obtained as 


tw 2 = -2 Pk 2 + fc(A;/p* - 1) 

• [A-jp* + 1 + P '- I ' i k- 1 ' 2 A* B (2 a;/ p - + 1 )~ 3/4 ]- 

(53) 

Stability in the region A* < —p*j 2 below the 
cellular boundary is determined by the real part of 
iu 2 . In that region, the principal value of the com- 
plex factor in Eq. (53) may be written in the form 
(A;/p* + 1)- 3/4 = [ - (A;/p* + l)]- 3/ V 3 -/ 4 , and 

thus the neutral stability condition lle{iu) 2 } = 0 
leads to an implicit equation for the (pulsating) neu- 
tral stability boundary A*(k,A* e] ,P). In terms of 
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the new pressure sensitivity parameter 6 defined by 
A* = -(p*/2)(l + 6), where 6 represents the neg- 
ative deviation, in units of p*/ 2 from the cellular 
boundary A* = -p*/ 2, this boundary is given by 

6 3/2 (3 + 6) -2 [(3 + 6)(l-6) + 8Pfc] 2 = a 1/2 /^> (54) 

w here oc = 4Aq ^ / p* - In the limit fc * oo, it is clear 
that that are two solutions of Eq. (54) given by 

6 = 0 (i.e., A* = -p*/2) and 6 1 + 2(1 + 2 Pk) 1/2 

(i.e., A'/p* (1 + 2Pk) 1/2 . Thus, the pulsating 

boundary is clearly multi-valued, as in the inviscid 
case (Fig. 3), with one branch approaching the cel- 
lular boundary and the other branch approaching 
the pulsating boundary for Aq = 0 (Fig. 2) in the 
limit of large fc. More generally, Eq. (54) may be 
rewritten as a cubic equation for the inverse relation 
fc(6) as 

64P 2 Jfc 3 + 16(3 + 6)(1 - 6)Pfc 2 

+ (3 + 6) 2 (1 - 6 ) 2 fc - a 1/2 ( 3 + 6 ) 2 6' 3/ 2 = 0 , 

(55) 

which is clearly seen to collapse to the previous in- 
viscid result 7 in the limit P -♦ 0. For arbitrary P, 
typical plots of fc(6) are shown in Figs. 4a— d, which, 
when rotated -90° so that the fc-axis is horizontal, 
is readily interpreted in the context of Figs. 1-3, 
where the lines A * = — p’/2 and A* = — p corre- 
spond to 6 = 0 and 6 = 1, respectively. It is clear 
that these curves asymptote to the lines 6 = 0 and 
6 = _i -f 2(1 + 2 Pfc) 1/2 as fc -* oo, where the lat- 
ter corresponds to the viscous pulsating boundary in 
the limit Aq — > 0. They cross the line 6 = 1, which 
corresponds to the inviscid pulsating boundary in 
the above limit, at fc 3 = a 1 / 2 /4P 2 . The fact that 
the pulsating boundary becomes C-shaped (in the 
rotated frame of reference) for A*q > 0 implies that 
steady, planar burning is intrinsically unstable for 
sufficiently small wavenumbers. In addition, since 
the portion within the C-shaped curve is the sta- 
ble region, any crossing of the C-shaped boundary 
from the stable to the unstable region corresponds 
to the onset of a pulsating instability. As Aq in- 
creases, the turning point of the C-shaped pulsating 
boundary shifts to larger values of fc. On the other 
hand, as Aq becomes small, the turning point shifts 
to small values of fc such that this point eventually 
leaves the 0(1) wavenumber region for which Eq. 
(54) are valid. Indeed, it turns out that the tran- 
sition to separated pulsating and cellular branches 
occurs as Aq decreases through 0(e l/2 ) values for 
intermediate 0(£^ 2 ) wavenumbers. 7 Thus, as Aq 
becomes small, the original pulsating and cellular 


boundaries are recovered in the 0(1) wavenumber 
regime, but as Aq becomes large, the original cel- 
lular boundary lies within the unstable region for 
0(1) wavenumbers and the basic solution becomes 
intrinsically unstable to oscillatory disturbances. 

Figures 4a-d. The pulsating hydrodynamic stability 
boundaries for k ~ 0(1) and A© ~ in the 

general viscous case. Figures are drawn for the val- 
ues p* = 1 and (a) P = .001; (b) P = .01; (c) 
P = .1; (d) P = 1.0. 


Pulsating Stability Boundary 

k As- ^'"X * - 0(1): p'-l.P- 0.001 



Pulsating Stability Boundary 

* A& - 0(€ W X * - 0(1 V p' - A P • 0.01 



9 



Pulsating Stability Boundary 

* A.-Olr'I'M-Omp- ul.r.Ol 



ter A*, we denote the two solution branches of Eq. 
(54), which correspond to the portions of Figure 4 
that lie to the left and to the right of the turning- 
point minimum, by (fc) and A"*- 0,1 (k), where 

the superscript **o” denotes, as before, the outer, 
or 0(1), wavenumber region and the superscripts 
“u” and “1” denote the upper and lower (rotate Fig- 
ure 4 by -90°) portions of the double-valued pul- 
sating boundary A* p {k). Along the upper branch, 
—p*/2 (i.e., b -» 0) as k -> oo, which can 
be matched with Eq. (36) since A* (/) -» -p * / 2 as 
kf — * 0. Similarly, A* < '°' 1 ^ — ♦ — p*( 1 + 2 Pfc) 1 ^ 2 (ie., 
& _1 + 2(1 4- 2Pk) 1/2 ) as fc -» oo, which clearly 

matches the viscous pulsating boundary given by the 
second of Eqs. (38) in the far outer wavenumber re- 
gion. As a result, a leading-order composite stabil- 
ity boundary spanning both the outer and far outer 
i wavenumber regions is given by 


Figure 4c 


Pulsating Stability Boundary 



Figure 4d 


Composite Neutral Stability Boundary 

A composite asymptotic solution for the neu- 
tral stability boundary in the regime A© ~ 0(e 1/4 ) 
is thus obtained by matching the cellular and pulsat- 
ing boundaries in the far outer wavenumber regime, 
where the former is given by Eq. (36) and the latter 
by the second of Eqs. (38), with the appropriate so- 
lution branch of Eq. (54) in the 0(1) wavenumber 
region. In particular, reverting back to the parame- 


f A; { °' u) (k)-p*/2 

I 2 p*p*P [1 + tfc(p* 7 + 2 p*P + 2 p*P)] 

+ 4p*P(l+p*Pefc) - [1— P( 6 fc)J(p *7 + 2/x*P) 

a; { 0 ' 1 ) Wi 


for A* % A* c , where A* p denotes the turning point 
calculated from Eq. (54) and the second term in the 
top expression has been expressed in terms of the 
outer wavenumber variable k. 

The composite stability boundary is shown in 
Fig. 5. Based on the above construction, the lower 
branch of Eq. (56) is a pulsating boundary for all 
wavenumbers, whereas the upper branch transitions 
from a pulsating boundary for 0(1) wavenumbers 
to a cellular boundary for 0(e _1 ) wavenumbers. In- 
deed, from Eq. (45), the size of the upper region 
of oscillatory instability, which is bounded below by 
the upper branch of the pulsating stability boundary 
and above by the region of nonoscillatory instability 

beyond the outer cellular boundary A* p* /2 for 

Aq = 0, shrinks to zero as k becomes large on the 
0(1) wavenumber scale. In this regime, the lack of a 
stable region for sufficiently small wavenumbers thus 
implies an intrinsic instability to long-wave pulsat- 


ing perturbations. 
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Composite Pulsating/ Cellular Stability Boundary 

A- @ =.5,p‘ =y = p* = l,P = .01.e = .005 





Figure 4. The composite pulsating/ cellular hydrody- 
namic stability boundary for Aq ~ 0(e 1,/4 ) in the 
general viscous case. 

Conclusion 

The present work further extends our recent for- 
mal treatment of hydrodynamic instability in liquid- 
propellant combustion. The analysis is based on a 
generalized Landau/Levich model in which the dy- 
namic motion of the liquid/gas interface, assumed 
to coincide with the combustion front, realistically 
possesses both a pressure and temperature sensi- 
tivity. In the present work, the fully viscous case 
was considered, thereby generalizing previous anal- 
yses in which either the viscosity of the fluid and/or 
the temperature sensitivity of the reaction rate was 
neglected. As in these preceding studies, the small- 
ness of the gas-to-liquid density ratio was used to de- 
fine a small parameter that allowed an asymptotic 
treatment of a rather complex dispersion relation. 
Specifically, it was again shown that in addition to 
the classical Landau, or cellular, stability bound- 
ary, there exists a pulsating hydrodynamic stability 
boundary as well. For sufficiently small values of the 
temperature-sensitivity parameter, there is a stable 
region between these two boundaries corresponding 
to a range of negative pressure sensitivities for which 
steady, planar burning is stable. 


As the pressure sensitivity decreases in mag- 
nitude, the cellular stability threshold is crossed, 
leading to classical Landau instability. Surface ten- 
sion, viscosity (both liquid and gas), and gravity 
are all stabilizing effects with respect to this type 
of instability. However, only gravity can stabilize 
small-wavenumber disturbances, and thus Landau 
instability becomes a long-wave instability in the 
reduced-gravity limit. Alternatively, as the pressure- 
sensitivity parameter increases in magnitude, the 
pulsating boundary is crossed, and liquid-propellant 
combustion becomes unstable to oscillatory pertur- 
bations. This type of hydrodynamic instability is 
more sensitive to the stabilizing effects of (liquid) 
viscosity than is the cellular boundary, but the stabi- 
lizing influence of viscosity does not extend to small 
wavenumber disturbances, and gravity turns out not 
to have a significant effect on this type of hydrody- 
namic instability. Consequently, for sufficiently large 
values of the temperature-sensitivity parameter, the 
pulsating boundary develops a turning point and 
becomes C-shaped. In this parameter regime, cor- 
responding to ratios of the temperature-to-pressure 
sensitivities of the order of 200 — 1000, steady, pla- 
nar combustion is intrinsically unstable to nonsteady 
long- wave perturbations. In that case, the pulsat- 
ing form of hydrodynamic instability is predicted to 
dominate, leading to large unsteady cells along the 
burning liquid/gas interface. 
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Nomenclature 


A 

Ap , A& 
hi 


e 

Ft 

9 

k 

n s 

P 

P, Pr 

9 

r 

t 

u 

V 

(x, y , z) 

7 

c 

< 

A 

P 

p 

<p3 

CJ 


burning rate 

pressure-, temperature-sensitivity 
coefficients 

coefficients in perturbation solution (i = 
1 , 2 , ... , 10 ) 
rate-of-strain tensor 
Froude number 

inverse Froude number (gravitational ac- 
celeration) 

perturbation wavenumber 

unit normal 

pressure 

Prandtl number 

quantity defined below Eq. (26) 

quantity defined below Eq. (26) 

time variable 

perturbation velocity vector 
velocity vector 
moving coordinate system 
surface-tension coefficient 
small bookkeeping parameter 
perturbation pressure 
gas-toliquid thermal diffusivity ratio 
gas-to-liquid viscosity ratio 
gas-to-liquid density ratio 
perturbation in location of gas/liquid inter- 
face 

location of gas/liquid interface 
complex perturbation frequency 


Subscripts, Superscripts: 

i inner wavenumber regime or integer vari- 
able 

/ far outer wavenumber regime 
l liquid 
9 gas 

o outer wavenumber regime 
* scaled quantity 
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